# This script will produce the tables in Avery & Brevoort, "The Subprime Crisis:  Is Government Housing Policy to Blame"

rm(list=ls())  # Clear the workspace

library(car)

# Create a log file
sink(file="Avery_Brevoort_Final_Estimations.log",
     append=FALSE,
     type="output",
     split=TRUE)

# This function rounds number to 3 digits and converts to a string
Trim <- function( x, digits=3 ) {
  outVal <- format( round(x,digits=digits) )
  return(outVal)
}

# This function converts p-values into stars for use in output table
Stars <- function( pValue ) {
  outval <- rep( "   ", times=length(pValue) )
  outval[ pValue<=0.1 ] <- "*  "
  outval[ pValue<=0.05 ] <- "** "
  outval[ pValue<=0.01 ] <- "***"
  return(outval)
}

# This function will be used to assemble regression output into output tables.
# j is the output of an OLS regression
# numPlaces provides the number of variables that will be included in the table (plus 1)
makeTable <- function( j, numPlaces ) {
  temp <- summary(j)$coefficients
  coefficients <- Trim( as.numeric(temp[2:numPlaces,1]), digits=3 )
  stars <- Stars( temp[2:numPlaces,4] )
  stderrs <- Trim( as.numeric(temp[2:numPlaces,2]), digits=3 )
  outVal <- paste( coefficients, stars, stderrs )
  return(outVal)
}


# Read in the tab-delimited data file
load("REStatSubmission_FRB.RData")

#######################################################################################
### Table 1 - Summary Statistics 
#######################################################################################

# Calculate the values in each row of the Summary Statistics table, Table 1
table1.output <- c("2008 Delinquency Rate", Trim(weighted.mean(myData$delrate,w=myData$COUNT), digits=2) )
table1.output <- rbind( table1.output, c("Center City", Trim(weighted.mean(myData$CENTER,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("Median House Value", Trim(weighted.mean(myData$CUML,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("Median Age", Trim(weighted.mean(myData$MEDAGE,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% 65 or Older", Trim(weighted.mean(myData$OVER65,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("Unemployment Rate", Trim(weighted.mean(myData$UNEMRT,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% Minority", Trim(weighted.mean(myData$minpct,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("Owner Occupancy Rate", Trim(weighted.mean(myData$ownpct,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("Vacancy Rate", Trim(weighted.mean(myData$vacpct,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2000 Relative Income", Trim(weighted.mean(myData$relinc,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("1990 Relative Income", Trim(weighted.mean(myData$relinc90,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Credit Score", Trim(weighted.mean(myData$SCORE2000,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004 Credit Score", Trim(weighted.mean(myData$SCORE2004,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% High PTI", Trim(weighted.mean(myData$hidti_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% High Rate", Trim(weighted.mean(myData$hirate_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% Piggyback", Trim(weighted.mean(myData$pig_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% No Income", Trim(weighted.mean(myData$zerinc_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% LMI", Trim(weighted.mean(myData$lowinc_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("% Below Median Income", Trim(weighted.mean(myData$midinc_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Dep Out Share", Trim(weighted.mean(myData$shr.dep.out_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Dep In Share", Trim(weighted.mean(myData$shr.dep.in_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Aff Out Share", Trim(weighted.mean(myData$shr.aff.out_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Aff In Share", Trim(weighted.mean(myData$shr.aff.in_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Credit Union Share", Trim(weighted.mean(myData$shr.cu_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Independent Share", Trim(weighted.mean(myData$shr.ind_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Dep Out Share", Trim(weighted.mean(myData$shr.dep.out_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Dep In Share", Trim(weighted.mean(myData$shr.dep.in_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Aff Out Share", Trim(weighted.mean(myData$shr.aff.out_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Aff In Share", Trim(weighted.mean(myData$shr.aff.in_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Credit Union Share", Trim(weighted.mean(myData$shr.cu_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Independent Share", Trim(weighted.mean(myData$shr.ind_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Dep Out Purchases", Trim(weighted.mean(myData$pur.dep.out_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Dep In Purchases", Trim(weighted.mean(myData$pur.dep.in_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Aff Out Purchases", Trim(weighted.mean(myData$pur.aff.out_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Aff In Purchases", Trim(weighted.mean(myData$pur.aff.in_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Credit Union Purchases", Trim(weighted.mean(myData$pur.cu_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 Independent Purchases", Trim(weighted.mean(myData$pur.ind_2001,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2001 GSE Sales", Trim(weighted.mean(myData$gse_2001,w=myData$COUNT),digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Dep Out Purchases", Trim(weighted.mean(myData$pur.dep.out_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Dep In Purchases", Trim(weighted.mean(myData$pur.dep.in_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Aff Out Purchases", Trim(weighted.mean(myData$pur.aff.out_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Aff In Purchases", Trim(weighted.mean(myData$pur.aff.in_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Credit Union Purchases", Trim(weighted.mean(myData$pur.cu_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 Independent Purchases", Trim(weighted.mean(myData$pur.ind_20046,w=myData$COUNT), digits=2) ) )
table1.output <- rbind( table1.output, c("2004-6 GSE Sales", Trim(weighted.mean(myData$gse_20046,w=myData$COUNT),digits=2) ) )

#######################################################################################
### Table 2 
#######################################################################################

# table2.a to table2.c estimate models using only MSA fixed effects.  
# These results will be used in calculating the within-MSA R-squared
table2.a <- lm(formula=delrate ~ as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & L11==1),
               model=FALSE)
table2.b <- lm(formula=hidti_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & L11==1),
               model=FALSE)
table2.c <- lm(formula=hirate_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & L11==1),
               model=FALSE)


# table2.1 to table2.5 estimate the models presented in Table 2
table2.1 <- lm(formula=delrate ~ shr.dep.out_2001 + shr.dep.in_2001 + shr.aff.out_2001 + shr.aff.in_2001 + shr.cu_2001 +
                                  pur.dep.out_2001 + pur.dep.in_2001 + pur.aff.out_2001 + pur.aff.in_2001 + pur.cu_2001 + pur.ind_2001 +
                                  gse_2001 + 
                                  SCORE2000 +
                                  relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                  as.factor(MSA),
                data=myData,
                weights=COUNT,
                subset=( MINOF3==1 & L11==1),
                model=TRUE)
table2.2 <- lm(formula=delrate ~ shr.dep.out_20046 + shr.dep.in_20046 + shr.aff.out_20046 + shr.aff.in_20046 + shr.cu_20046 +
                                 pur.dep.out_20046 + pur.dep.in_20046 + pur.aff.out_20046 + pur.aff.in_20046 + pur.cu_20046 + pur.ind_20046 +
                                 gse_20046 +
                                 SCORE2004 +
                                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & L11==1),
               model=TRUE)
table2.3 <- lm(formula=delrate ~ shr.dep.out_20046 + shr.dep.in_20046 + shr.aff.out_20046 + shr.aff.in_20046 + shr.cu_20046 +
                                 pur.dep.out_20046 + pur.dep.in_20046 + pur.aff.out_20046 + pur.aff.in_20046 + pur.cu_20046 + pur.ind_20046 +
                                 gse_20046 + 
                                 SCORE2004 +
                                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                 hidti_20046 + hirate_20046 + pig_20046 + lowinc_20046 + midinc_20046 + zerinc_20046 +
                                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & L11==1),
               model=TRUE)
table2.4 <- lm(formula=hidti_20046 ~ shr.dep.out_2001 + shr.dep.in_2001 + shr.aff.out_2001 + shr.aff.in_2001 + shr.cu_2001 +
                                     pur.dep.out_2001 + pur.dep.in_2001 + pur.aff.out_2001 + pur.aff.in_2001 + pur.cu_2001 + pur.ind_2001 +
                                     gse_2001 + 
                                     SCORE2000 +
                                     relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                     as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & L11==1),
               model=TRUE)
table2.5 <- lm(formula=hirate_20046 ~ shr.dep.out_2001 + shr.dep.in_2001 + shr.aff.out_2001 + shr.aff.in_2001 + shr.cu_2001 +
                                      pur.dep.out_2001 + pur.dep.in_2001 + pur.aff.out_2001 + pur.aff.in_2001 + pur.cu_2001 + pur.ind_2001 +
                                      gse_2001 + 
                                      SCORE2000 +
                                      relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                      as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & L11==1),
               model=TRUE)

# Assemble the coefficients into an output table
table2.1.coefs <- makeTable( table2.1, numPlaces=13 )
table2.2.coefs <- makeTable( table2.2, numPlaces=13 )
table2.3.coefs <- makeTable( table2.3, numPlaces=13 )
table2.4.coefs <- makeTable( table2.4, numPlaces=13 )
table2.5.coefs <- makeTable( table2.5, numPlaces=13 )
table2.names <- dimnames( summary(table2.1)$coefficients )[[1]][2:13]
table2.coefs <- cbind( table2.names, table2.1.coefs, table2.2.coefs, table2.3.coefs, table2.4.coefs, table2.5.coefs )

# This next block of code calculates the within-MSA R-squared for each estimation, based
# on the SSEs from each estimation and from the models estimated using only the MSA-level
# fixed effects.
table2.a.sse <- weighted.mean(x=table2.a$residuals^2, w=table2.a$weights)*length(table2.a$weights)
table2.b.sse <- weighted.mean(x=table2.b$residuals^2, w=table2.b$weights)*length(table2.b$weights)
table2.c.sse <- weighted.mean(x=table2.c$residuals^2, w=table2.c$weights)*length(table2.c$weights)
table2.1.sse <- weighted.mean(x=table2.1$residuals^2, w=table2.1$weights)*length(table2.1$residuals)
table2.2.sse <- weighted.mean(x=table2.2$residuals^2, w=table2.2$weights)*length(table2.2$residuals)
table2.3.sse <- weighted.mean(x=table2.3$residuals^2, w=table2.3$weights)*length(table2.3$residuals)
table2.4.sse <- weighted.mean(x=table2.4$residuals^2, w=table2.4$weights)*length(table2.4$residuals)
table2.5.sse <- weighted.mean(x=table2.5$residuals^2, w=table2.5$weights)*length(table2.5$residuals)
table2.1.wR2 <- Trim( 1 - table2.1.sse/table2.a.sse, digits=3 )
table2.2.wR2 <- Trim( 1 - table2.2.sse/table2.a.sse, digits=3 )
table2.3.wR2 <- Trim( 1 - table2.3.sse/table2.a.sse, digits=3 )
table2.4.wR2 <- Trim( 1 - table2.4.sse/table2.b.sse, digits=3 )
table2.5.wR2 <- Trim( 1 - table2.5.sse/table2.c.sse, digits=3 )
table2.wR2 <- cbind( "R-squared", table2.1.wR2, table2.2.wR2, table2.3.wR2, table2.4.wR2, table2.5.wR2 )

# Report the number of observations in each estimation
table2.1.obs <- length(table2.1$model[,1])
table2.2.obs <- length(table2.2$model[,2])
table2.3.obs <- length(table2.3$model[,2])
table2.4.obs <- length(table2.4$model[,2])
table2.5.obs <- length(table2.5$model[,2])
table2.obs <- cbind( "Observations", table2.1.obs, table2.2.obs, table2.3.obs, table2.4.obs, table2.5.obs)

# Report the mean of the dependent variable in each estimation
table2.1.ybar <- Trim( weighted.mean(x=table2.1$model[,1], w=table2.1$weights), digits=3 )
table2.2.ybar <- Trim( weighted.mean(x=table2.2$model[,1], w=table2.2$weights), digits=3 )
table2.3.ybar <- Trim( weighted.mean(x=table2.3$model[,1], w=table2.3$weights), digits=3 )
table2.4.ybar <- Trim( weighted.mean(x=table2.4$model[,1], w=table2.4$weights), digits=3 )
table2.5.ybar <- Trim( weighted.mean(x=table2.5$model[,1], w=table2.5$weights), digits=3 )
table2.ybar <- cbind( "Mean y", table2.1.ybar, table2.2.ybar, table2.3.ybar, table2.4.ybar, table2.5.ybar)

# Conduct the three statistical tests (t1, t2, and t3) for each specification
table2.1.t1 <- linearHypothesis(table2.1,"shr.dep.out_2001 = shr.dep.in_2001",test="F")
table2.1.t2 <- linearHypothesis(table2.1,"shr.aff.out_2001 = shr.aff.in_2001",test="F")
table2.1.t3 <- linearHypothesis(table2.1,c("shr.dep.out_2001 = shr.dep.in_2001","shr.aff.out_2001 = shr.aff.in_2001"),test="F")
table2.2.t1 <- linearHypothesis(table2.2,"shr.dep.out_20046 = shr.dep.in_20046",test="F")
table2.2.t2 <- linearHypothesis(table2.2,"shr.aff.out_20046 = shr.aff.in_20046",test="F")
table2.2.t3 <- linearHypothesis(table2.2,c("shr.dep.out_20046 = shr.dep.in_20046","shr.aff.out_20046 = shr.aff.in_20046"),test="F")
table2.3.t1 <- linearHypothesis(table2.3,"shr.dep.out_20046 = shr.dep.in_20046",test="F")
table2.3.t2 <- linearHypothesis(table2.3,"shr.aff.out_20046 = shr.aff.in_20046",test="F")
table2.3.t3 <- linearHypothesis(table2.3,c("shr.dep.out_20046 = shr.dep.in_20046","shr.aff.out_20046 = shr.aff.in_20046"),test="F")
table2.4.t1 <- linearHypothesis(table2.4,"shr.dep.out_2001 = shr.dep.in_2001",test="F")
table2.4.t2 <- linearHypothesis(table2.4,"shr.aff.out_2001 = shr.aff.in_2001",test="F")
table2.4.t3 <- linearHypothesis(table2.4,c("shr.dep.out_2001 = shr.dep.in_2001","shr.aff.out_2001 = shr.aff.in_2001"),test="F")
table2.5.t1 <- linearHypothesis(table2.5,"shr.dep.out_2001 = shr.dep.in_2001",test="F")
table2.5.t2 <- linearHypothesis(table2.5,"shr.aff.out_2001 = shr.aff.in_2001",test="F")
table2.5.t3 <- linearHypothesis(table2.5,c("shr.dep.out_2001 = shr.dep.in_2001","shr.aff.out_2001 = shr.aff.in_2001"),test="F")

# Format the output of the statistical tests
table2.1.t1.out <- paste( Trim(table2.1.t1$F[2],digits=2), Stars(table2.1.t1[6][2,1]) )
table2.2.t1.out <- paste( Trim(table2.2.t1$F[2],digits=2), Stars(table2.2.t1[6][2,1]) )
table2.3.t1.out <- paste( Trim(table2.3.t1$F[2],digits=2), Stars(table2.3.t1[6][2,1]) )
table2.4.t1.out <- paste( Trim(table2.4.t1$F[2],digits=2), Stars(table2.4.t1[6][2,1]) )
table2.5.t1.out <- paste( Trim(table2.5.t1$F[2],digits=2), Stars(table2.5.t1[6][2,1]) )
table2.1.t2.out <- paste( Trim(table2.1.t2$F[2],digits=2), Stars(table2.1.t2[6][2,1]) )
table2.2.t2.out <- paste( Trim(table2.2.t2$F[2],digits=2), Stars(table2.2.t2[6][2,1]) )
table2.3.t2.out <- paste( Trim(table2.3.t2$F[2],digits=2), Stars(table2.3.t2[6][2,1]) )
table2.4.t2.out <- paste( Trim(table2.4.t2$F[2],digits=2), Stars(table2.4.t2[6][2,1]) )
table2.5.t2.out <- paste( Trim(table2.5.t2$F[2],digits=2), Stars(table2.5.t2[6][2,1]) )
table2.1.t3.out <- paste( Trim(table2.1.t3$F[2],digits=2), Stars(table2.1.t3[6][2,1]) )
table2.2.t3.out <- paste( Trim(table2.2.t3$F[2],digits=2), Stars(table2.2.t3[6][2,1]) )
table2.3.t3.out <- paste( Trim(table2.3.t3$F[2],digits=2), Stars(table2.3.t3[6][2,1]) )
table2.4.t3.out <- paste( Trim(table2.4.t3$F[2],digits=2), Stars(table2.4.t3[6][2,1]) )
table2.5.t3.out <- paste( Trim(table2.5.t3$F[2],digits=2), Stars(table2.5.t3[6][2,1]) )
table2.t1 <- cbind( "T1", table2.1.t1.out, table2.2.t1.out, table2.3.t1.out, table2.4.t1.out, table2.5.t1.out)
table2.t2 <- cbind( "T2", table2.1.t2.out, table2.2.t2.out, table2.3.t2.out, table2.4.t2.out, table2.5.t2.out)
table2.t3 <- cbind( "T3", table2.1.t3.out, table2.2.t3.out, table2.3.t3.out, table2.4.t3.out, table2.5.t3.out)

# Assemble the pieces of the output table
# This variable will contain all of the information displayed in table 2
table2.output <- rbind( table2.coefs, table2.obs, table2.wR2, table2.ybar, table2.t1, table2.t2, table2.t3 )

#######################################################################################
### Table 3 - Threshold Estimations for CRA
#######################################################################################

# table3.a to table3.e estimate models using only MSA fixed effects.  
# These results will be used in calculating the within-MSA R-squared
table3.a <- lm(formula=delrate ~ as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=FALSE)
table3.b <- lm(formula=hidti_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=FALSE)
table3.c <- lm(formula=hirate_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=FALSE)
table3.d <- lm(formula=shr.dep.in_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=FALSE)
table3.e <- lm(formula=pur.dep.in_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=FALSE)

# table3.1 to table3.6 estimate the models presented in Table 3
table3.1 <- lm(formula=delrate ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                 SCORE2000 +
                                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)
table3.2 <- lm(formula=delrate ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                 SCORE2004 +
                                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                 hidti_20046 + hirate_20046 + pig_20046 + lowinc_20046 + midinc_20046 + zerinc_20046 +
                                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)
table3.3 <- lm(formula=hidti_20046 ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                     SCORE2000 +
                                     UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                     as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)
table3.4 <- lm(formula=hirate_20046 ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                      SCORE2000 +
                                      UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                      as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)
table3.5 <- lm(formula=shr.dep.in_20046 ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                          SCORE2000 +
                                          UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                          as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)
table3.6 <- lm(formula=pur.dep.in_20046 ~ cra_75 + cra_76 + cra_77 + cra_78 + cra_79 + cra_80 + cra_81 + cra_82 + cra_83 +
                                          SCORE2000 +
                                          UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                                          as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=75 & relinc<85),
               model=TRUE)

# Assemble the coefficients into an output table
table3.1.coefs <- makeTable( table3.1, numPlaces=10 )
table3.2.coefs <- makeTable( table3.2, numPlaces=10 )
table3.3.coefs <- makeTable( table3.3, numPlaces=10 )
table3.4.coefs <- makeTable( table3.4, numPlaces=10 )
table3.5.coefs <- makeTable( table3.5, numPlaces=10 )
table3.6.coefs <- makeTable( table3.6, numPlaces=10 )
table3.names <- dimnames( summary(table3.1)$coefficients )[[1]][2:10]
table3.coefs <- cbind( table3.names, table3.1.coefs, table3.2.coefs, table3.3.coefs, table3.4.coefs, table3.5.coefs, table3.6.coefs )

# This next block of code calculates the within-MSA R-squared for each estimation, based
# on the SSEs from each estimation and from the models estimated using only the MSA-level
# fixed effects.
table3.a.sse <- weighted.mean(x=table3.a$residuals^2, w=table3.a$weights)*length(table3.a$weights)
table3.b.sse <- weighted.mean(x=table3.b$residuals^2, w=table3.b$weights)*length(table3.b$weights)
table3.c.sse <- weighted.mean(x=table3.c$residuals^2, w=table3.c$weights)*length(table3.c$weights)
table3.d.sse <- weighted.mean(x=table3.d$residuals^2, w=table3.d$weights)*length(table3.d$weights)
table3.e.sse <- weighted.mean(x=table3.e$residuals^2, w=table3.e$weights)*length(table3.e$weights)
table3.1.sse <- weighted.mean(x=table3.1$residuals^2, w=table3.1$weights)*length(table3.1$residuals)
table3.2.sse <- weighted.mean(x=table3.2$residuals^2, w=table3.2$weights)*length(table3.2$residuals)
table3.3.sse <- weighted.mean(x=table3.3$residuals^2, w=table3.3$weights)*length(table3.3$residuals)
table3.4.sse <- weighted.mean(x=table3.4$residuals^2, w=table3.4$weights)*length(table3.4$residuals)
table3.5.sse <- weighted.mean(x=table3.5$residuals^2, w=table3.5$weights)*length(table3.5$residuals)
table3.6.sse <- weighted.mean(x=table3.6$residuals^2, w=table3.6$weights)*length(table3.6$residuals)
table3.1.wR2 <- Trim( 1 - table3.1.sse/table3.a.sse, digits=3 )
table3.2.wR2 <- Trim( 1 - table3.2.sse/table3.a.sse, digits=3 )
table3.3.wR2 <- Trim( 1 - table3.3.sse/table3.b.sse, digits=3 )
table3.4.wR2 <- Trim( 1 - table3.4.sse/table3.c.sse, digits=3 )
table3.5.wR2 <- Trim( 1 - table3.5.sse/table3.d.sse, digits=3 )
table3.6.wR2 <- Trim( 1 - table3.6.sse/table3.e.sse, digits=3 )
table3.wR2 <- cbind( "R-squared", table3.1.wR2, table3.2.wR2, table3.3.wR2, table3.4.wR2, table3.5.wR2, table3.6.wR2 )

# Report the number of observations in each estimation
table3.1.obs <- length(table3.1$model[,1])
table3.2.obs <- length(table3.2$model[,2])
table3.3.obs <- length(table3.3$model[,2])
table3.4.obs <- length(table3.4$model[,2])
table3.5.obs <- length(table3.5$model[,2])
table3.6.obs <- length(table3.6$model[,2])
table3.obs <- cbind( "Observations", table3.1.obs, table3.2.obs, table3.3.obs, table3.4.obs, table3.5.obs, table3.6.obs)

# Report the mean of the dependent variable in each estimation
table3.1.ybar <- Trim( weighted.mean(x=table3.1$model[,1], w=table3.1$weights), digits=3 )
table3.2.ybar <- Trim( weighted.mean(x=table3.2$model[,1], w=table3.2$weights), digits=3 )
table3.3.ybar <- Trim( weighted.mean(x=table3.3$model[,1], w=table3.3$weights), digits=3 )
table3.4.ybar <- Trim( weighted.mean(x=table3.4$model[,1], w=table3.4$weights), digits=3 )
table3.5.ybar <- Trim( weighted.mean(x=table3.5$model[,1], w=table3.5$weights), digits=3 )
table3.6.ybar <- Trim( weighted.mean(x=table3.6$model[,1], w=table3.6$weights), digits=3 )
table3.ybar <- cbind( "Mean Y", table3.1.ybar, table3.2.ybar, table3.3.ybar, table3.4.ybar, table3.5.ybar, table3.6.ybar)

# Conduct the three statistical tests (t1 and t2) for each specification
table3.1.t1 <- linearHypothesis(table3.1,"2*cra_79 = cra_78 + cra_80",test="F")
table3.1.t2 <- linearHypothesis(table3.1,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")
table3.2.t1 <- linearHypothesis(table3.2,"2*cra_79 = cra_78 + cra_80",test="F")
table3.2.t2 <- linearHypothesis(table3.2,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")
table3.3.t1 <- linearHypothesis(table3.3,"2*cra_79 = cra_78 + cra_80",test="F")
table3.3.t2 <- linearHypothesis(table3.3,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")
table3.4.t1 <- linearHypothesis(table3.4,"2*cra_79 = cra_78 + cra_80",test="F")
table3.4.t2 <- linearHypothesis(table3.4,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")
table3.5.t1 <- linearHypothesis(table3.5,"2*cra_79 = cra_78 + cra_80",test="F")
table3.5.t2 <- linearHypothesis(table3.5,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")
table3.6.t1 <- linearHypothesis(table3.6,"2*cra_79 = cra_78 + cra_80",test="F")
table3.6.t2 <- linearHypothesis(table3.6,"8*cra_79 = cra_75 + cra_76 + cra_77 + cra_78 + cra_80 + cra_81 + cra_82 + cra_83")

# Format the output of the statistical tests
table3.1.t1.out <- paste( Trim(table3.1.t1$F[2],digits=2), Stars(table3.1.t1[6][2,1]) )
table3.2.t1.out <- paste( Trim(table3.2.t1$F[2],digits=2), Stars(table3.2.t1[6][2,1]) )
table3.3.t1.out <- paste( Trim(table3.3.t1$F[2],digits=2), Stars(table3.3.t1[6][2,1]) )
table3.4.t1.out <- paste( Trim(table3.4.t1$F[2],digits=2), Stars(table3.4.t1[6][2,1]) )
table3.5.t1.out <- paste( Trim(table3.5.t1$F[2],digits=2), Stars(table3.5.t1[6][2,1]) )
table3.6.t1.out <- paste( Trim(table3.6.t1$F[2],digits=2), Stars(table3.6.t1[6][2,1]) )
table3.1.t2.out <- paste( Trim(table3.1.t2$F[2],digits=2), Stars(table3.1.t2[6][2,1]) )
table3.2.t2.out <- paste( Trim(table3.2.t2$F[2],digits=2), Stars(table3.2.t2[6][2,1]) )
table3.3.t2.out <- paste( Trim(table3.3.t2$F[2],digits=2), Stars(table3.3.t2[6][2,1]) )
table3.4.t2.out <- paste( Trim(table3.4.t2$F[2],digits=2), Stars(table3.4.t2[6][2,1]) )
table3.5.t2.out <- paste( Trim(table3.5.t2$F[2],digits=2), Stars(table3.5.t2[6][2,1]) )
table3.6.t2.out <- paste( Trim(table3.6.t2$F[2],digits=2), Stars(table3.6.t2[6][2,1]) )
table3.t1 <- cbind( "T1", table3.1.t1.out, table3.2.t1.out, table3.3.t1.out, table3.4.t1.out, table3.5.t1.out, table3.6.t1.out)
table3.t2 <- cbind( "T2", table3.1.t2.out, table3.2.t2.out, table3.3.t2.out, table3.4.t2.out, table3.5.t2.out, table3.6.t2.out)

# Assemble the pieces of the output table
# This variable will contain all of the information displayed in table 2
table3.output <- rbind( table3.coefs, table3.obs, table3.wR2, table3.ybar, table3.t1, table3.t2 )

#######################################################################################
### Table 4 - Threshold Estimations for GSE Income
#######################################################################################

# table4.a to table4.d estimate models using only MSA fixed effects.  
# These results will be used in calculating the within-MSA R-squared
table4.a <- lm(formula=delrate ~ as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=FALSE)
table4.b <- lm(formula=hidti_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=FALSE)
table4.c <- lm(formula=hirate_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=FALSE)
table4.d <- lm(formula=gse_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=FALSE)

# table4.1 to table4.5 estimate the models presented in Table 4
table4.1 <- lm(formula=delrate ~ gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_89 + gsei_90 + gsei_91 + gsei_92 + gsei_93 +
                 SCORE2000 +
                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=TRUE)
table4.2 <- lm(formula=delrate ~ gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_89 + gsei_90 + gsei_91 + gsei_92 + gsei_93 +
                 SCORE2004 +
                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                 hidti_20046 + hirate_20046 + pig_20046 + lowinc_20046 + midinc_20046 + zerinc_20046 +
                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=TRUE)
table4.3 <- lm(formula=hidti_20046 ~ gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_89 + gsei_90 + gsei_91 + gsei_92 + gsei_93 +
                 SCORE2000 +
                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=TRUE)
table4.4 <- lm(formula=hirate_20046 ~ gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_89 + gsei_90 + gsei_91 + gsei_92 + gsei_93 +
                 SCORE2000 +
                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=TRUE)
table4.5 <- lm(formula=gse_20046 ~ gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_89 + gsei_90 + gsei_91 + gsei_92 + gsei_93 +
                 SCORE2000 +
                 UNEMRT + MEDAGE + OVER65 + CENTER + minpct + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=85 & relinc<95 & minpct<=30),
               model=TRUE)

# Assemble the coefficients into an output table
table4.1.coefs <- makeTable( table4.1, numPlaces=10 )
table4.2.coefs <- makeTable( table4.2, numPlaces=10 )
table4.3.coefs <- makeTable( table4.3, numPlaces=10 )
table4.4.coefs <- makeTable( table4.4, numPlaces=10 )
table4.5.coefs <- makeTable( table4.5, numPlaces=10 )
table4.names <- dimnames( summary(table4.1)$coefficients )[[1]][2:10]
table4.coefs <- cbind( table4.names, table4.1.coefs, table4.2.coefs, table4.3.coefs, table4.4.coefs, table4.5.coefs )

# This next block of code calculates the within-MSA R-squared for each estimation, based
# on the SSEs from each estimation and from the models estimated using only the MSA-level
# fixed effects.
table4.a.sse <- weighted.mean(x=table4.a$residuals^2, w=table4.a$weights)*length(table4.a$weights)
table4.b.sse <- weighted.mean(x=table4.b$residuals^2, w=table4.b$weights)*length(table4.b$weights)
table4.c.sse <- weighted.mean(x=table4.c$residuals^2, w=table4.c$weights)*length(table4.c$weights)
table4.d.sse <- weighted.mean(x=table4.d$residuals^2, w=table4.d$weights)*length(table4.d$weights)
table4.1.sse <- weighted.mean(x=table4.1$residuals^2, w=table4.1$weights)*length(table4.1$residuals)
table4.2.sse <- weighted.mean(x=table4.2$residuals^2, w=table4.2$weights)*length(table4.2$residuals)
table4.3.sse <- weighted.mean(x=table4.3$residuals^2, w=table4.3$weights)*length(table4.3$residuals)
table4.4.sse <- weighted.mean(x=table4.4$residuals^2, w=table4.4$weights)*length(table4.4$residuals)
table4.5.sse <- weighted.mean(x=table4.5$residuals^2, w=table4.5$weights)*length(table4.5$residuals)
table4.1.wR2 <- Trim( 1 - table4.1.sse/table4.a.sse, digits=3 )
table4.2.wR2 <- Trim( 1 - table4.2.sse/table4.a.sse, digits=3 )
table4.3.wR2 <- Trim( 1 - table4.3.sse/table4.b.sse, digits=3 )
table4.4.wR2 <- Trim( 1 - table4.4.sse/table4.c.sse, digits=3 )
table4.5.wR2 <- Trim( 1 - table4.5.sse/table4.d.sse, digits=3 )
table4.wR2 <- cbind( "R-squared", table4.1.wR2, table4.2.wR2, table4.3.wR2, table4.4.wR2, table4.5.wR2 )

# Report the number of observations in each estimation
table4.1.obs <- length(table4.1$model[,1])
table4.2.obs <- length(table4.2$model[,2])
table4.3.obs <- length(table4.3$model[,2])
table4.4.obs <- length(table4.4$model[,2])
table4.5.obs <- length(table4.5$model[,2])
table4.obs <- cbind( "Observations", table4.1.obs, table4.2.obs, table4.3.obs, table4.4.obs, table4.5.obs)

# Report the mean of the dependent variable in each estimation
table4.1.ybar <- Trim( weighted.mean(x=table4.1$model[,1], w=table4.1$weights), digits=3 )
table4.2.ybar <- Trim( weighted.mean(x=table4.2$model[,1], w=table4.2$weights), digits=3 )
table4.3.ybar <- Trim( weighted.mean(x=table4.3$model[,1], w=table4.3$weights), digits=3 )
table4.4.ybar <- Trim( weighted.mean(x=table4.4$model[,1], w=table4.4$weights), digits=3 )
table4.5.ybar <- Trim( weighted.mean(x=table4.5$model[,1], w=table4.5$weights), digits=3 )
table4.ybar <- cbind( "Mean Y", table4.1.ybar, table4.2.ybar, table4.3.ybar, table4.4.ybar, table4.5.ybar)

# Conduct the three statistical tests (t1 and t2) for each specification
table4.1.t1 <- linearHypothesis(table4.1,"2*gsei_89 = gsei_88 + gsei_90",test="F")
table4.1.t2 <- linearHypothesis(table4.1,"8*gsei_89 = gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_90 + gsei_91 + gsei_92 + gsei_93")
table4.2.t1 <- linearHypothesis(table4.2,"2*gsei_89 = gsei_88 + gsei_90",test="F")
table4.2.t2 <- linearHypothesis(table4.2,"8*gsei_89 = gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_90 + gsei_91 + gsei_92 + gsei_93")
table4.3.t1 <- linearHypothesis(table4.3,"2*gsei_89 = gsei_88 + gsei_90",test="F")
table4.3.t2 <- linearHypothesis(table4.3,"8*gsei_89 = gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_90 + gsei_91 + gsei_92 + gsei_93")
table4.4.t1 <- linearHypothesis(table4.4,"2*gsei_89 = gsei_88 + gsei_90",test="F")
table4.4.t2 <- linearHypothesis(table4.4,"8*gsei_89 = gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_90 + gsei_91 + gsei_92 + gsei_93")
table4.5.t1 <- linearHypothesis(table4.5,"2*gsei_89 = gsei_88 + gsei_90",test="F")
table4.5.t2 <- linearHypothesis(table4.5,"8*gsei_89 = gsei_85 + gsei_86 + gsei_87 + gsei_88 + gsei_90 + gsei_91 + gsei_92 + gsei_93")

# Format the output of the statistical tests
table4.1.t1.out <- paste( Trim(table4.1.t1$F[2],digits=2), Stars(table4.1.t1[6][2,1]) )
table4.2.t1.out <- paste( Trim(table4.2.t1$F[2],digits=2), Stars(table4.2.t1[6][2,1]) )
table4.3.t1.out <- paste( Trim(table4.3.t1$F[2],digits=2), Stars(table4.3.t1[6][2,1]) )
table4.4.t1.out <- paste( Trim(table4.4.t1$F[2],digits=2), Stars(table4.4.t1[6][2,1]) )
table4.5.t1.out <- paste( Trim(table4.5.t1$F[2],digits=2), Stars(table4.5.t1[6][2,1]) )
table4.1.t2.out <- paste( Trim(table4.1.t2$F[2],digits=2), Stars(table4.1.t2[6][2,1]) )
table4.2.t2.out <- paste( Trim(table4.2.t2$F[2],digits=2), Stars(table4.2.t2[6][2,1]) )
table4.3.t2.out <- paste( Trim(table4.3.t2$F[2],digits=2), Stars(table4.3.t2[6][2,1]) )
table4.4.t2.out <- paste( Trim(table4.4.t2$F[2],digits=2), Stars(table4.4.t2[6][2,1]) )
table4.5.t2.out <- paste( Trim(table4.5.t2$F[2],digits=2), Stars(table4.5.t2[6][2,1]) )
table4.t1 <- cbind( "T1", table4.1.t1.out, table4.2.t1.out, table4.3.t1.out, table4.4.t1.out, table4.5.t1.out)
table4.t2 <- cbind( "T2", table4.1.t2.out, table4.2.t2.out, table4.3.t2.out, table4.4.t2.out, table4.5.t2.out)

# Assemble the pieces of the output table
# This variable will contain all of the information displayed in table 2
table4.output <- rbind( table4.coefs, table4.obs, table4.wR2, table4.ybar, table4.t1, table4.t2 )

#######################################################################################
### Table 5 - Threshold Estimations for GSE Minority
#######################################################################################

# table5.a to table5.d estimate models using only MSA fixed effects.  
# These results will be used in calculating the within-MSA R-squared
table5.a <- lm(formula=delrate ~ as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=FALSE)
table5.b <- lm(formula=hidti_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=FALSE)
table5.c <- lm(formula=hirate_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=FALSE)
table5.d <- lm(formula=gse_20046 ~ as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=FALSE)

# table5.1 to table5.5 estimate the models presented in Table 5
table5.1 <- lm(formula=delrate ~ gsem_34 + gsem_33 + gsem_32 + gsem_31 + gsem_30 + gsem_29 + gsem_28 + gsem_27 + gsem_26 +
                 SCORE2000 +
                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=TRUE)
table5.2 <- lm(formula=delrate ~ gsem_34 + gsem_33 + gsem_32 + gsem_31 + gsem_30 + gsem_29 + gsem_28 + gsem_27 + gsem_26 +
                 SCORE2004 +
                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + CUML + ownpct + vacpct +
                 hidti_20046 + hirate_20046 + pig_20046 + lowinc_20046 + midinc_20046 + zerinc_20046 +
                 as.factor(MSA),
               data=myData,
               weights=COUNT,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=TRUE)
table5.3 <- lm(formula=hidti_20046 ~ gsem_34 + gsem_33 + gsem_32 + gsem_31 + gsem_30 + gsem_29 + gsem_28 + gsem_27 + gsem_26 +
                 SCORE2000 +
                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=TRUE)
table5.4 <- lm(formula=hirate_20046 ~ gsem_34 + gsem_33 + gsem_32 + gsem_31 + gsem_30 + gsem_29 + gsem_28 + gsem_27 + gsem_26 +
                 SCORE2000 +
                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=TRUE)
table5.5 <- lm(formula=gse_20046 ~ gsem_34 + gsem_33 + gsem_32 + gsem_31 + gsem_30 + gsem_29 + gsem_28 + gsem_27 + gsem_26 +
                 SCORE2000 +
                 relinc + relinc90 + UNEMRT + MEDAGE + OVER65 + CENTER + CUML + ownpct + vacpct +
                 as.factor(MSA),
               data=myData,
               weights=loans_20046,
               subset=( MINOF3==1 & relinc>=90 & relinc<120 & minpct>=25 & minpct<35),
               model=TRUE)

# Assemble the coefficients into an output table
table5.1.coefs <- makeTable( table5.1, numPlaces=10 )
table5.2.coefs <- makeTable( table5.2, numPlaces=10 )
table5.3.coefs <- makeTable( table5.3, numPlaces=10 )
table5.4.coefs <- makeTable( table5.4, numPlaces=10 )
table5.5.coefs <- makeTable( table5.5, numPlaces=10 )
table5.names <- dimnames( summary(table5.1)$coefficients )[[1]][2:10]
table5.coefs <- cbind( table5.names, table5.1.coefs, table5.2.coefs, table5.3.coefs, table5.4.coefs, table5.5.coefs )

# This next block of code calculates the within-MSA R-squared for each estimation, based
# on the SSEs from each estimation and from the models estimated using only the MSA-level
# fixed effects.
table5.a.sse <- weighted.mean(x=table5.a$residuals^2, w=table5.a$weights)*length(table5.a$weights)
table5.b.sse <- weighted.mean(x=table5.b$residuals^2, w=table5.b$weights)*length(table5.b$weights)
table5.c.sse <- weighted.mean(x=table5.c$residuals^2, w=table5.c$weights)*length(table5.c$weights)
table5.d.sse <- weighted.mean(x=table5.d$residuals^2, w=table5.d$weights)*length(table5.d$weights)
table5.1.sse <- weighted.mean(x=table5.1$residuals^2, w=table5.1$weights)*length(table5.1$residuals)
table5.2.sse <- weighted.mean(x=table5.2$residuals^2, w=table5.2$weights)*length(table5.2$residuals)
table5.3.sse <- weighted.mean(x=table5.3$residuals^2, w=table5.3$weights)*length(table5.3$residuals)
table5.4.sse <- weighted.mean(x=table5.4$residuals^2, w=table5.4$weights)*length(table5.4$residuals)
table5.5.sse <- weighted.mean(x=table5.5$residuals^2, w=table5.5$weights)*length(table5.5$residuals)
table5.1.wR2 <- Trim( 1 - table5.1.sse/table5.a.sse, digits=3 )
table5.2.wR2 <- Trim( 1 - table5.2.sse/table5.a.sse, digits=3 )
table5.3.wR2 <- Trim( 1 - table5.3.sse/table5.b.sse, digits=3 )
table5.4.wR2 <- Trim( 1 - table5.4.sse/table5.c.sse, digits=3 )
table5.5.wR2 <- Trim( 1 - table5.5.sse/table5.d.sse, digits=3 )
table5.wR2 <- cbind( "R-squared", table5.1.wR2, table5.2.wR2, table5.3.wR2, table5.4.wR2, table5.5.wR2 )

# Report the number of observations in each estimation
table5.1.obs <- length(table5.1$model[,1])
table5.2.obs <- length(table5.2$model[,2])
table5.3.obs <- length(table5.3$model[,2])
table5.4.obs <- length(table5.4$model[,2])
table5.5.obs <- length(table5.5$model[,2])
table5.obs <- cbind( "Observations", table5.1.obs, table5.2.obs, table5.3.obs, table5.4.obs, table5.5.obs)

# Report the mean of the dependent variable in each estimation
table5.1.ybar <- Trim( weighted.mean(x=table5.1$model[,1], w=table5.1$weights), digits=3 )
table5.2.ybar <- Trim( weighted.mean(x=table5.2$model[,1], w=table5.2$weights), digits=3 )
table5.3.ybar <- Trim( weighted.mean(x=table5.3$model[,1], w=table5.3$weights), digits=3 )
table5.4.ybar <- Trim( weighted.mean(x=table5.4$model[,1], w=table5.4$weights), digits=3 )
table5.5.ybar <- Trim( weighted.mean(x=table5.5$model[,1], w=table5.5$weights), digits=3 )
table5.ybar <- cbind( "Mean Y", table5.1.ybar, table5.2.ybar, table5.3.ybar, table5.4.ybar, table5.5.ybar)

# Conduct the three statistical tests (t1 and t2) for each specification
table5.1.t1 <- linearHypothesis(table5.1,"2*gsem_29 = gsem_28 + gsem_30",test="F")
table5.1.t2 <- linearHypothesis(table5.1,"8*gsem_29 = gsem_26 + gsem_27 + gsem_28 + gsem_30 + gsem_31 + gsem_32 + gsem_33 + gsem_34")
table5.2.t1 <- linearHypothesis(table5.2,"2*gsem_29 = gsem_28 + gsem_30",test="F")
table5.2.t2 <- linearHypothesis(table5.2,"8*gsem_29 = gsem_26 + gsem_27 + gsem_28 + gsem_30 + gsem_31 + gsem_32 + gsem_33 + gsem_34")
table5.3.t1 <- linearHypothesis(table5.3,"2*gsem_29 = gsem_28 + gsem_30",test="F")
table5.3.t2 <- linearHypothesis(table5.3,"8*gsem_29 = gsem_26 + gsem_27 + gsem_28 + gsem_30 + gsem_31 + gsem_32 + gsem_33 + gsem_34")
table5.4.t1 <- linearHypothesis(table5.4,"2*gsem_29 = gsem_28 + gsem_30",test="F")
table5.4.t2 <- linearHypothesis(table5.4,"8*gsem_29 = gsem_26 + gsem_27 + gsem_28 + gsem_30 + gsem_31 + gsem_32 + gsem_33 + gsem_34")
table5.5.t1 <- linearHypothesis(table5.5,"2*gsem_29 = gsem_28 + gsem_30",test="F")
table5.5.t2 <- linearHypothesis(table5.5,"8*gsem_29 = gsem_26 + gsem_27 + gsem_28 + gsem_30 + gsem_31 + gsem_32 + gsem_33 + gsem_34")

# Format the output of the statistical tests
table5.1.t1.out <- paste( Trim(table5.1.t1$F[2],digits=2), Stars(table5.1.t1[6][2,1]) )
table5.2.t1.out <- paste( Trim(table5.2.t1$F[2],digits=2), Stars(table5.2.t1[6][2,1]) )
table5.3.t1.out <- paste( Trim(table5.3.t1$F[2],digits=2), Stars(table5.3.t1[6][2,1]) )
table5.4.t1.out <- paste( Trim(table5.4.t1$F[2],digits=2), Stars(table5.4.t1[6][2,1]) )
table5.5.t1.out <- paste( Trim(table5.5.t1$F[2],digits=2), Stars(table5.5.t1[6][2,1]) )
table5.1.t2.out <- paste( Trim(table5.1.t2$F[2],digits=2), Stars(table5.1.t2[6][2,1]) )
table5.2.t2.out <- paste( Trim(table5.2.t2$F[2],digits=2), Stars(table5.2.t2[6][2,1]) )
table5.3.t2.out <- paste( Trim(table5.3.t2$F[2],digits=2), Stars(table5.3.t2[6][2,1]) )
table5.4.t2.out <- paste( Trim(table5.4.t2$F[2],digits=2), Stars(table5.4.t2[6][2,1]) )
table5.5.t2.out <- paste( Trim(table5.5.t2$F[2],digits=2), Stars(table5.5.t2[6][2,1]) )
table5.t1 <- cbind( "T1", table5.1.t1.out, table5.2.t1.out, table5.3.t1.out, table5.4.t1.out, table5.5.t1.out)
table5.t2 <- cbind( "T2", table5.1.t2.out, table5.2.t2.out, table5.3.t2.out, table5.4.t2.out, table5.5.t2.out)

# Assemble the pieces of the output table
# This variable will contain all of the information displayed in table 2
table5.output <- rbind( table5.coefs, table5.obs, table5.wR2, table5.ybar, table5.t1, table5.t2 )

# print out the tables that appear in the paper
print(table1.output)
print(table2.output)
print(table3.output)
print(table4.output)
print(table5.output)

sink() # Close log file

